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• Abstract. We revisit the issue of relaxation to thermal equilibrium in the so- 

C/2 . called "sheet model", i.e., particles in one dimension interacting by attractive forces 

independent of their separation. We show that this relaxation may be very clearly 
detected and characterized by following the evolution of order parameters defined by 
appropriately normalized moments of the phase space distribution which probe its 
' O ' entanglement in space and velocity coordinates. For a class of quasi-stationary states 

which result from the violent relaxation of rectangular waterbag initial conditions, 
^ I characterized by their virial ratio Rq, we show that relaxation occurs on a time scale 

which (i) scales approximately linearly in the particle number N , and (ii) shows also 
, a strong dependence on i?0: with quasi-stationary states from colder initial conditions 

^ ' relaxing much more rapidly. The temporal evolution of the order parameter may be 

■ well described by a stretched exponential function. We study finally the correlation of 

, the relaxation times with the amplitude of fluctuations in the relaxing quasi-stationary 

' states, as well as the relation between temporal and ensemble averages. 
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1. Introduction 

The so-called "sheet model" is an interesting toy model for the study of self-gravitating 
systems, or more generally of systems with long-range interactions. It is simply the one 
dimensional (ID) generalisation of Newtonian gravity, consisting of particles interacting 
by attractive forces independent of their separation (or, equivalently, infinite parallel 
planes embedded in three dimensions interacting via Newtonian gravity). Because the 
particle trajectories are exactly integrable between crossings, it has the nice feature that 
its numerical integration can be performed with an accuracy limited only by machine 
precision. It has been the subject of (mostly numerical) study in the literature for 
several decades (see, e.g., [1] for a review of the literature on the model) following 
earher analytical studies [21 13] . A fundamental question about this system — and more 
generally for any system with long-range interactions — is whether they relax to the 
statistical equilibrium calculated in the microcanonical or canonical ensemble. For this 
model the latter were first calculated exactly, for any particle number A^, by Rybicki [1]. 
The literature on this model — which we will discuss in greater detail in our conclusions 
section below — is marked by differing results (or, rather, interpretation of results) from 
different groups, and even some controversy. Work by two groups in the eighties (see, e.g. 
[5] for a summary) led to the conclusion that relaxation could not be observed, except 
perhaps in some special cases. Studies by two other groups over a decade ago P [7] 
found results indicating relaxation, and [6] gave a determination of the dependence 
of the characteristic time. However doubts about the interpretation of these latter 
results as establishing relaxation to equilibrium were raised by further analysis [H |9] . 
In more recent work [lOl [1] clear evidence for relaxation in a version of the model in 
which there are different particle masses has been found, but the dependence on A^ has 
not been determinec||| The mechanism of relaxation (if it indeed takes place) in these 
models remains, as in other long-range interacting systems, very poorly understood, 
and a basic subject of research in the statistical mechanics of long-range interacting 
systems (for recent reviews see e.g. [HI [l5]). In this article we report an essentially 
numerical study of relaxation in the single mass sheet model. We introduce a simple 
but, as we will see, very useful tool for the characterisation of the long-time evolution 
and relaxation of the system. This tool allows us to resolve some outstanding issues 
about the relaxation in this system, and, in particular, to establish more definitively 
both that relaxation does indeed occur and the scaling with particle number of the time 
characterizing it. We consider a broader range of initial conditions, which allows us to 
establish also dependences of relaxation on these. We also study the fluctuations — both 
in time and over realizations of the initial conditions — about the average macroscopic 
evolution of the system, showing phenomenologically the correlation of their amplitude 
with the lifetimes of the intermediate "quasi-stationary" states. 

We will discuss in greater detail in our conclusions the relation of our results to 
those in the previous literature, but it is useful at the outset to say a little more about 

I Other variants of the model have also been studied in [HI [121 US] • 
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the more general context of this study. In recent years there has been considerable 
interest in the statistical mechanics of long range interactions (see, e.g., [161 [13 IT^). 
stimulated by the need to understand the physics of various laboratory systems with 
interactions of this kind, as well as by the more classical case of self-gravitating matter 
relevant in astrophysics and cosmology. In this context one toy model in particular, and 
various variants of it, has been much studied: the Hamiltonian Mean Field (HMF) 
model (see e.g. fl8\ [T9| 120] and references therein), which is a model of particles 
on a circle interacting by a cosine potential. Its study has shown that it shares 
many of the qualitative features well documented in the most studied of realistic long- 
range interacting systems — self-gravitating systems in astrophysics — and believed 
to be generic in such systems. Starting from generic initial conditions, the system 
evolves rapidly (by "violent relaxation" ) to a virialized macroscopically stationary state. 
These states — commonly referred to in the more recent literature as "quasi-stationary 
states" (QSS) — are out of equilibrium states, which are described theoretically in the 
framework of Vlasov equation (more usually referred to as the "collisionless Boltzmann 
equations" in the astrophysical literature). On much longer time scales an evolution 
towards the true thermal equilibrium (i.e. that determined by the maximization of the 
Boltzmann entropy in a mean field approximation) is postulated. For realistic systems 
— such as Newtonian gravity in three dimensions — - it is very difficult numerically to 
simulate the evolution on sufficiently long time scales to probe the relaxation. Studies 
in the literature (see e.g. [211 1221 [231 [211 [25]) provide some results but give still a very 
limited characterization and understanding of it. The HMF model has the particular 
feature that the potential energy of any particle may be expressed as a function of its 
(angular) position and the mean potential energy due to all particles — it is for this 
reason that it is "mean-field" — so that the calculation of the forces in a system with 
particles requires only of order operations (rather than N"^ in a typical long-range 
interacting system). Further the force is continuous at zero separation, so that the 
difficulties associated in the case of gravity with the regulation of the potential at small 
scales are avoided. This allows the regime of relaxation to be accessed numerically even 
for quite large particle numbers. The study of [18] found a scaling of the relaxation time 
in proportion to iV^''' (but see also e.g. [I9] which finds indications of longer lifetimes 
for other initial conditions). 

It can clearly be of interest to study different toy models, to determine in particular 
features which are indeed generic. The "sheet model" is probably the oldest toy model 
of long range interactions — it was first explored in astrophysics as a toy model for self- 
gravitating systems in three dimensions — and is also, arguably, closer to reality than 
the HMF which is constrained on a circle. It has, further, as mentioned above the nice 
feature that it numerical integration can be performed up to machine precision. Despite 
this, the results concerning its dynamics and relaxation are less clearly determined than 
for the HMF, and the literature on the subject has, as we have discussed above, been 
marked by some controversy and results showing that the model has, apparently, some 
very peculiar behaviours — rapid relaxation to equilibrium for some classes of initial 
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states EZ], persistent phase space structures impeding relaxation to QSS [2H1 12^ . 
macroscopically chaotic behaviour in the long time evolution |30] — which indicate 
that it might not be a very useful toy model (in that its behaviours are perhaps non- 
generic). In this article our main conclusion is that, 1) by using appropriate diagnostics 
of the macroscopic evolution and 2) by extending simulations to sufficiently large 
and/or averaging over sufficiently large numbers of realization, one finds behaviour in 
this toy model very similar in crucial respects to that in the HMF: to a very good 
first approximation a generic initial configuration relaxes to a long-lived QSS, and then 
relaxes to its statistical equilibrium at sufficiently long time. This latter phase can be 
characterized apparently by a single time-scale, with the evolution of the order parameter 
during relaxation well fit by a simple function (in our case a better fit is obtained using 
a simple stretched exponential, rather than a hyperbolic tangent in the HMF as in [T8]). 
On the other hand the dependence of this time scale, linearly proportional to the 
number of particles A^, is different to that found in [18] for the HMF. This latter result, 
however, applies to spatially homogeneous states which, in the HMF, can occur due to 
the periodicity of the system. Relaxation which is slower than linear in N is expected 
in this shown using analysis of kinetic equations (see, e.g. contributions of P.H. 

Chavanis, and of F. Bouchet and J. Barre in [IT]). 

The article is organised as follows. In the next section we recall the basic definitions 
of the model, and relevant results on its statistical equilibrium. We then explain 
the choice of the macroscopic parameters ("order parameters") we choose to monitor 
the evolution of the system. In the following section we first describe our numerical 
simulations and the initial conditions we study, and then give our results. In presenting 
them we give first results for single realizations, and then use temporal averages and 
finally ensemble averages to derive the scaling with of the relaxation time. This is 
followed by further study of the fiuctuations about the average behaviours of the order 
parameters. Considering both temporal fiuctuations and those in the ensemble, which 
we show to be very consistent with one another, we observe the correlation between 
their amplitude in the QSS and the observed relaxation time. In the conclusion sections 
we return to a more detailed discussion of the previous literature, presenting further 
results which allow one to understand the reasons for the divergence in conclusions in 
certain cases. 

2. The sheet model 

We first recall the model and fix our notation. We next summarize the results of [^ on 
statistical equilibrium, and then explain the rationale for our choice of "order parameter" 
in our study. 
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2.1. Definitions 

We consider identical (equal mass) point particles in one dimension interacting by an 
attractive force independent of separation, i.e., the force f{x) on a particle at coordinate 
position X exerted by a particle at the origin is given by 



where g is the coupling. Equivalently it is the pair interaction derived from the pair 
potential 



which satisfies the ID Poisson equation for a point source, = 2gSD{x) (where Sd is 
the Dirac delta function). Comparing with the three dimensional (3D) Poisson equation 
shows the equivalence with the case of an infinitely thin plane of infinite extent and 
surface mass density S = g/2'n'G, which explains the widely used name "sheet model". 
We will work in the one dimensional language, referring to "particles". For a particle 
at coordinate position x on the real axis the total force F{x) acting on it is thus 



where Ny{x) and N^{x) are, respectively, the number of particles with coordinates 
greater than or less than x (i.e. the force on a given particle is proportional to the 
difference in the number of particles on its right and its left). 

To specify fully the dynamics we must prescribe what happens when two particles 
arrive at the same point. Since the force is bounded as the separation goes to zero, 
the natural physical prescription for the ID model is that the particles simply cross (i.e. 
pass through one another). In one dimension, however, this is equivalent, up to a change 
in particle labels, to a hard elastic collision, as such a collision (of equal mass particles) 
simply results in an exchange of their velocities. Thus, up to particle labels, the sheet 
model for equal masses is equivalent to one in which particles experience always the same 
spatially constant force Eq. ([3]) and simply exchange velocities when they "collide" . As 
has been noted in some previous studies of models of this kind [3T] it is convenient to 
exploit this equivalence in numerical simulation, as will be described below. 

In contrast to Newtonian gravity in three dimensions, the pair potential ([2]) is 
positive and diverges at large separations, so that particles cannot escape from the 
system to infinity. It has therefore no particular interest to enclose the system in a finite 
box, and indeed such a confinement is not necessary in order to define the statistical 
equilibrium (in contrast to three dimensions). We will consider therefore always open 
boundary conditions. Likewise the fact that the potential has no divergence at short 
distances means that there is no equivalent of the so-called "gravo-thermal collapse" in 
three dimensions. 




(1) 




(2) 




(3) 
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2.2. Thermal equilibrium 

It has been shown by Rybicki [4J that the statistical equihbrium for this model can be 
derived exactly in the microcanonanical ensemble, for any A^. We will study here the 
N ^ oo limit (at fixed energy and mass, see jl] for full derivation), in which case the 
phase space distribution function (i.e. mass per unit phase space volume) becomes: 

/e.(x,.) = ^-^e ^sech- (4) 

where a and A are the characteristic scales of velocity and length, and M is the total 
mass of the system. It is straightforward to verify that 

= ^ (6) 

3M ^ ^ 

^ - '■'^ 

where E is the total energy of the system, which allows one to calculate feq{x,v) 
explicitly as a function of M and E (and g) only. 

As is typical of long-range systems, the statistical equilibrium is thus characterised 
by a space independent Maxwellian velocity distribution and an inhomogeneous spatial 
distribution. The same solution is recovered in the canonical ensemble. Thus, differently 
to many long-range systems (including gravity in three dimensions) the two ensembles 
are completely equivalent. This behaviour is associated also with the absence of 
microcanonical phase transitions which may arise in such systems. 

This equilibrium solution in the continuum limit can be most easily derived by 
simply maximizing the Boltzmann entropy at fixed mass and energy, using the mean- 
field expression for the energy: 

E=-J v'^ f{x,v)dxdv + - J f{x,v)ip{\x — x'\)f{x' ,v')dxdx'dvdv' . (7) 

This procedure gives simply 

^ 1/2 2 

Ux, v) = Ae-^[^(-)+^l = !^e"^p(x) (8) 
where (l){x) is the mean field potential 

(j){x)= I i){\x-x'\)f{x\v')dx'dv' . (9) 



and p{x) the associated mass density profile, which is therefore the solution to 

^[Inp(x)] = -2gPp{x). (10) 

It is simple to verify that the expression Eq. (jlj) results, with appropriate identification 
of constants and choice of units. 
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2. 3. Order parameters for relaxation 

To monitor relaxation to equilibrium it is possible in principle to simply study the full 
distribution function as a function of time. In practice relaxation is extremely slow (in 
the characteristic time units of the system) and only accessible numerically for relatively 
small numbers (of order a thousand) particles, which makes the comparison of the full 
function subtle (because of finite fluctuations). In the previous literature various 
methods have been used — statistical tests on the velocity and spatial distributions 
(e.g. [261 E]) cind analyses based on the evolution of particle energies coarse-grained 
in time [6]). Here instead we study the evolution primarily using appropriately chosen 
macroscopic parameters, i.e., "order parameters" which take, in general, distinct values 
in and out of equilibrium. This is somewhat analogous to the approach used in the 
study of the HMF, where the magnetisation of the system plays the role of an order 
parameter used to characterize the evolution out of equilibrium (see e.g. [18J). Once the 
expected behaviour of the macroscopic parameter is identified, a more detailed analysis 
involving the distribution functions can be used to confirm that the system has indeed 
fully relaxed. We will see that, with the choice of parameter we make, it turns out that 
the single macroscopic parameter is sometimes a better indicator of relaxation that the 
full density or velocity distributions, and that indeed some of the controversial results 
in the previous literature may easily be sorted out using the tools used here. 

An evident property of the distribution function Eq. (jl]) is that it is separable in its 
spatial and velocity coordinates. It is simple to show, as we will now verify, that it is in 
fact also the unique stationary solution of the Vlasov equation which is separable. Thus 
if the system is, during its evolution, very close at any time to a stationary solution 
of the Vlasov equation (which describes the collisionless limit) any parameter probing 
the degree of separabihty of the distribution function would be expected to be a useful 
indicative measure. This leads us to consider order parameters which are simply suitably 
normalized moments of the distribution function. 

Let us first verify the result on separability. The Vlasov equation for the model is 



which can be conveniently expressed in terms of the mean field potential 0(x) satisfying 
the Poisson equation with p{x) as source , i.e.. 



with an appropriate boundary condition from Eq. (fT2!) [e.g. a{x — )■ +oo) = gM, where 
M is the total mass]. Seeking a solution which is both stationary and separable we take 





where a{x) is the mean field acceleration, i.e., 

a{x) = g I sgn{x' — x)f{x' ^v)dx'dv 




(12) 




(13) 




(14) 
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and thus obtain, on substitution, 

«(.„).. „.^ = _p(,).„W.^, (15) 

Given that a{x) ^ everywhere (except at the single point which divides the mass in 
two) we can write this as 1§| 

1 dejv) _ 1 dpjx) ^ ^^g^ 

V ■ 6{v) dv a{x) ■ p{x) dx 

for any region of x where p{x) ^ 0. It follows that both sides are equal to a constant, 
C say, and therefore 

B(y) = Q^e-'i'^\ (17) 

The right hand side gives 

Differentiating with respect to x and using the Poisson equation for the mean field 
Eq. (IT^ . this gives exactly the same equation used above to determine the equilibrium 
solution for p{x). The expression Eq. (jl]) is indeed therefore the only stationary separable 
solution of the Vlasov equation. 

Given this observation we define the following order parameters: 

1 (19) 





1 


\v\ 








1^) 



for non-zero a and /3, where 
1 ^ 

= (20) 

i=\ 

and Ui is the value of the parameter u for the z-th particle. By construction these 
quantities are zero in thermal equilibrium. While a finite number of such moments 
can of course be zero in a QSS with a non-separable distribution function, we expect 
them generically to be non-zero in such states. We will use here both 0ii and 022- As 
detailed in the next section, we will consider both their temporal evolution in single 
realizations of our initial conditions, as well as averages of these temporal evolutions. 
These averages will be performed in two different ways: by averaging over a finite 
temporal window in a single realization, and by averaging over independent realizations 
of the initial conditions. Further we will consider the evolution as a function of time of 
the fluctuations of and 022 with respect to these averages. 

§ Note that this is not true in the HMF model, as the magnetization (which determines the acceleration) 
can indeed be zero everywhere when the stationary state is spatially uniform. This is a result of the 
periodic nature of the system. In this case there may thus exist QSS which are separable, uniform in 
space but with a non- maxwellian velocity distribution. Such QSS are indeed observed and have been 
extensively studied in this model (see e.g. [18]). 
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3. Numerical simulations 



3.1. Algorithm 

As remarked above in Sec. I2.H it is convenient for the numerical integration of the model 
to exchange particles' labels when they cross, which is equivalent to treating them as if 
they undergo an elastic collision in which they exchange their velocities when they meet. 
The force on each particle is then constant in space and time [and given by Eq. ([3])], 
and the numerical algorithm must simply determine, at any time, the next crossing 
which occurs, and then exchange the velocities of the "colliding" particles. The optimal 
way to treat this kind of problem is, as has been pointed out and discussed in detail in 
[32] . by using a so-called "heap-based" algorithm, which uses an object called a "heap" 
to store in an ordered way the next crossing times of the pairs (see [32] for details). 
This algorithm requires a number of operations of order log(A^) to determine the next 
crossing (rather than of order N for the evident direct algorithm in which one calculates 
and compares directly at each step the next crossing of each of the — 1 pairs) . Given 
that the number of crossings per particle per unit time grows in proportion to A^, the 
simulation time thus grows in proportion to N'^log{N). 

Because the particle trajectories are integrated exactly, the only limit on the 
accuracy of the numerical integration is thus the numerical precision. As is common 
practice we will use the total energy (which is conserved in the continuum model) as 
a control parameter. For the longest simulations we report the error in total energy of 
the order of 10"^%. 



3.2. Initial conditions 

We will consider principally a simple class of spatially uniform initial conditions (IC), 
generated by randomly distributing the A^ particles on a finite interval. As initial velocity 
distribution we will consider both the case that initial velocities are zero ("cold IC") 
and the case that this distribution is also uniform in a finite interval. The latter are thus 
random samplings of a particular class of "waterbag" initial conditions in phase space 
(i.e. in which the phase space density in equal in the region in which it is non-zero), 
while the cold case can be considered as the limit in which the width of the velocity 
distribution goes to zero. In Fig.[T]the phase space distribution for a typical IC is shown. 

These IC may be characterized solely by the particle number A^ and a single 
parameter characterizing the waterbag. Rather than the width of the velocity 
distribution or phase space density, it is convenient to choose the parameter 
characterizing the waterbag to be the dimensionless initial virial ratio: 

Ro^'^ (21) 

where Tq is the initial average total kinetic energy and Uq is the initial average total 
potential energy. By "average" here we mean that the values of Tq and Uq are those 
calculated for the theoretical waterbag configuration. When we consider, as we do 
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Figure 1. A rectangular waterbag initial condition in phase space sampled with 
N = 400 particles, for a initial virial Rq — 1 (see text for definition of units). 



below, different realizations of the initial conditions at fixed and -Rq there are of 
course finite N fluctuations about these values of Tq and Uq. The latter correspond 
then to their average values in the ensemble of initial conditions at fixed N and Rq. We 
note that this ensemble of initial conditions is clearly not a subset of the microcanonical 
ensemble because there are finite N fluctuations also about the average total energy. 
As these fluctuations are small, however, we will assume below that the evolution of 
such an ensemble of initial conditions represents well that of an ensemble of such initial 
conditions with exactly the average energjjjj]. 

We remark on a particularity of the cold IC which we will return to below. In the 
limit — )■ oo, the evolution from this IC becomes singular at a finite time: an element 
of mass initially at coordinate position xq feels a force —2gpoXo, where po is the initial 
mass density; all particles are in free-fall under a force proportional to their distance, 
and therefore arrive at the origin at the same time, producing a density singularity. For 
a finite system, the corresponding behaviour is associated to the existence of a periodic 
oscillating mode when the particles are initially equally spaced (i.e. on a regular lattice). 
This "breathing mode" of such a cold system has been discussed in [33]. While there is 
no such mode in a three dimensional system of a finite number of particles, the N ^ oo 
limit of a cold spherical initial condition has the analogous singularity (see [S^ for a 
detailed discussion). 

II It is simple to show, given that the particle positions and velocities are both randomly sampled from 
a PDF which is uniform in a finite interval, that the relative fluctuations in Uq and Tq scale as 
for large N . An exact calculation shows, for example, that at iV = 100, the normalized variance of C/q, 
which corresponds to that in the energy for the case i?o = 0, is « 0.05. This means that the typical 
amplitude of the fluctuation in the energy for cold IC at = 100 is of order the difference between the 
mean energy of cold IC and IC with Rq — 0.1. 
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3.3. Units and coordinates 

For convenience we choose our coordinate system such that the centre of the mass of the 
system hes at x = and is at rest (i.e. after distributing the particles as described we 
add a spatial translation and constant velocity to all particles so that these conditions 
are satisfied). 

We make the following choice of units: we set the particle mass m and the coupling 
g to unity, and take L = N. This corresponds to a mass (and particle) density of unity, 
and a time unit equal to the dynamical time: 



which is the characteristic time for the system's evolution under the mean field forces 
(the mean field forces, of order Ng, moves a system particle of mass m over the system 
size L on this time-scale), tdyn also coincides with the time of the singularity noted 
above in the smooth limit of the cold IC. 

4. Results 

The difficulty in this study of relaxation, as in such a study for any long-range system, 
is that one is interested in studying large N systems — so that finite deviations 
from the mean-field behaviour are small — on a time scale which grows rapidly with 
(typically, one expects, in proportion to some power of N). Because of numerical 
limitations, particularly strong because of the computational cost of integrating a long- 
range interaction, it is in practice often difficult to arrive at definitive conclusions. In 
the case of gravity in three dimensions, notably, numerical studies exist (see references 
above) but they give only a very incomplete characterization and understanding of 
relaxation. As we have discussed in the introduction one of the attractive features of 
the HMF model is that, because of its mean field nature, the numerical cost of the force 
calculation is of order A^, allowing much larger particle numbers — A^ ~ 10^ — 10^ [18] 
— to be simulated on the relevant long time scales than is feasible in other cases. The 
principal reason why the early literature on the sheet model was marked by controversy 
on the question of relaxation is simply, as we will discuss further below, that such 
relaxation could not be observed on the required time scales for systems sufficiently large 
so that the finite A^ fluctuations were sufficiently small to allow the clear identification 
of the average behaviours. The study of [6], taking advantage of the greater numerical 
resources available already in the nineties, detected relaxation for A^ ~ 10^ from specific 
waterbag configurations and found a scaling of the relaxation time, over a small range 
in A^, linear in A^. This result was obtained, however, by doing a time average of 
their chosen diagnostic over a very broad time window (of order 10^, only an order 
of magnitude less than the typical relaxation time for the cases explored), and its 
solidity has been placed in question in subsequent work [8l|9]. Exploiting the increase in 
numerical power since then, and aided greatly by the diagnostics we have defined in the 




(22) 
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previous section, we report here [^results showing relaxation for systems with ~ 10^. 
Further we obtain our results for the scaling of the relaxation time by doing ensemble 
averages (over realizations of the initial conditions) without time averages. 

4.I. Temporal evolution of order parameters 
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Figure 2. Evolution as a function of time of the virial ratio (left panel), (j>ii (middle 
panel) and 4)22 (right panel), in a simulation with N = 100 particles and initial virial 
ratio Ro = 0- We observe that the system virializes on a time scale of a few tens of 
dynamical times; further the behaviours of 0ii and (/>22 indicate, as typically expected 
in long-range interacting systems, a subsequent evolution in a long-lived QSS which 
eventually relaxes to thermal equilibrium (in which 01 1 = = (f>22 = on average). 



Shown in Fig. |2] is the evolution of the virial ratio R, and the order parameters 
011 and 022 in a single realization of a system with N = 100, and Rq = 0. Note that 
the time axis (as it will be invariably here) is logarithmic. In Fig. [3] are plotted the 
same quantities for = 400. While the fluctuations are very large, particularly in the 
first case, one can make out that there are, as expected, two stages in the macroscopic 
evolution probed by these parameters: a first stage (t < 100) of "violent relaxation" 
during which all quantities (and notably the virial ratio) fiuctuates strongly before 
settling down to behaviours which appear to fiuctuate about a well defined average, 
and specifically about unity for the virial ratio. The averages of the parameters (pu 
and 022 are clearly non-zero on a much longer time scale than that characterizing the 
virialization. These non-zero averages, which appear to be approximately the same in 
each case for the two different A^, appear to remain roughly stable until at least about 
10^ — 10^, after which both 0ii and 022 start to evolve towards zero. The time scale at 
which the evolution sets in is clearly significantly shorter for = 100. This behaviour 
should indicate, as we have discussed above, the relaxation to statistical equilibrium. 

These behaviours can be seen more clearly by averaging in a temporal window, 
of width small compared to the characteristic times scales of this apparent evolution. 

% Evolution of = 10^ particles to t = 10^ requires about 20 minutes on a single processor; thus, 
given the scaling with iV^ logiV of the computational cost per unit physical time, and a linear growth 
(see below) in the relaxation time itself, simulation times of order several weeks are required for the 
most rapidly relaxing case with N = 10"^. Our largest N results are ensemble averages over systems 
with N = 800, obtained by running simultaneously on a large number of work stations. 
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Figure 3. Evolution as a function of time of the same three parameters as in Fig. 
[21 in a simulation with N = 400 and initial virial ratio Rq = 0. We see the same 
qualitative behaviours as in the previous figure, except that fluctuations are of lower 
amplitude and the QSS phase appears to persist longer. 




Figure 4. Evolution as a function of time of the same three quantities and the same 
initial conditions as in Fig. [5J but with these quantities now averaged in a time 
window of width At = 10 as described in text. The averaging makes the interpretation 
given in Fig. [2] much clearer: once virialized the system stays in a long-lived QSS and 
eventually relaxes to thermal equilibrium. 



Shown in Fig. H] and |5] are the same quantities for the same simulations, but now each 
point represents the average over one hundred time shces, equally spaced in a window 
of width At = 10 centred on the given time. 




Figure 5. Evolution as a function of time of the same quantities and for the same 
initial conditions with N = 400 as in Fig. [31 but with these quantities now averaged 
in a time window of width At = 10. Comparing to the previous figure (same quantities 
for A'^ = 100) we see clearly that the QSS persists for longer. 



These behaviours are thus clearly indicative of the evolution expected, which is that 
believed to be typical of long-range interacting systems: violent relaxation brings one on 
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a short time scale to a QSS, as a result of a mean field dynamics described by the coupled 
Vlasov-Poisson equations (and this independent of A^). On longer, N — dependent times 
scales, one relaxes to the mean-field equilibrium, given in this case by Eq. (jl]). That 
the decay to zero of 0ii and 022 does indeed correspond to relaxation to the statistical 
equilibrium of Eq. (jlj) can be tested in further detail. Fig. [6] shows the velocity and 
space distributions for Rq = and = 400 particles, averaged again over a time 
window of width At = 10, at t = 10^ and t = 10^. The continuous lines correspond to 
Eq. (H]), clearly in very good agreement at the later time, and very different in the QSS 
phase. We have also checked (but do not show here) the agreement of the distribution 
of particle energies. These results indicate that (pu and 022 are very good diagnostics 
of the evolution towards equilibrium: indeed below we will see that they are typically 
more discriminating of relaxation than the full density and velocity distributions. 
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Figure 6. Density profiles (top, left and right panels) and velocity distributions 
(bottom, left and right panels) aX t — 10"^ (left) and t = 10^ (right), in a simulation 
with N — 400 particles started from initial conditions with i?o = 0. The quantities 
are averaged in a temporal window of width At exactly as in the previous two figures. 
The continuous lines are the expected distributions at thermal equilibrium, Eq. ([4]). 



4-2. Dependence on initial virial ratio 

Shown in Fig. [7] are the evolution of 0ii for = 100 and N = 400 starting now from 
four different values of Rq, as indicated (0,0.1,0.5,1) The results are averaged again 
in a time window of width At = 10. Note that results for N = 100, which extend up 
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to t = 10®, indicate that the evolution towards the statistical equilibrium is really a 
relaxation to a definitive equilibrium behaviour, i.e., which is stable and persists. This 
is further confirmed by Fig. [8] which shows the spatial and velocity distributions for 
the case Rq = and N = 100 at t = 10® (with the same time averaging window as used 
above in Fig. |6]). Fig. [7]showv that there are, nevertheless, very significant fiuctuations 
in and 022- These could indicate significant macroscopic, but stochastic, deviations 
from the equilibrium persisting over very significant times (see [30]). We will present 
evidence below that they are, as one would expect, finite effects, with an amplitude 
which decreases as N increases. 
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Figure 7. Evolution as a function of time of 0ii, averaged in a time window of width 
A ~ 10, for simulations with N — 100 particles (left panel) and N — 400 particles 
(right panel), starting from initial conditions with the indicated values of Rq. The 
time scale for relaxation of a QSS clearly depends not just on N but on the details of 
this state. 
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Figure 8. Density profile (left panel) and velocity distributions (right panel), averaged 
in a time window of width A = 10 centered at t = 10^, for a simulation with N = 100 
particles started from a virial ratio Rq = 0. The continuous lines are the expected 
distributions at thermal equilibrium, Eq. (j4|). 



We observe in Fig. [7|that, as expected, the QSS resulting from violent relaxation 
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is clearly different depending on the initial condition, with very different values of 
011. Further the relaxation toward equilibrium is evident in most cases, but at a time 
which depends not only on A^, but also on Rq (or the intermediate QSS state). More 
specifically, the smaller is Rq the shorter is the lifetime. Indeed for i?o = 1 we just see 
the onset of the relaxation for the case N = 100, but do not see it at all for = 400. 
For N = 100 there is a difference of a factor of about one hundred in the time at which 
relaxation appears to becomes clearly visible in the cases Rq = 1 and Rq = 0. In 
the respect we remark that earlier studies have not considered this kind of cold initial 
condition, in which relaxation occurs more rapidly. 



4-3. Estimation of N dependence using ensemble average 

Let us focus now on the N dependence of the relaxation. We wish to determine the 
scaling with of the characteristic time for relaxation, at a fixed value of the initial virial 
ratio. Given the very significant noise in the order parameters at the particle numbers 
we can simulate numerically up to times on which relaxation occurs to do so we must 
average out these fluctuations. This can be done using either a time average on a single 
realization (as above) or an average over realizations (or possibly some combination of 
both). 
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Figure 9. Evolution as a function of time of 0ii for N ~ 100 particles (left panel) 
and N = 800 particles (right panel) in simulations starting from three different 
realizations of initial conditions with i?o = 0. An average in a temporal window 
of width At = 10 has been used in all cases. Despite the time averaging there are still 
significant fluctuations which limit the precision of the determination of the time scale 
for relaxation. 



Shown in Fig. |9] is a plot of the evolution of 0ii in three different realizations for 
= 100 and A^ = 800 and Rq = 0, up to t = 10^. The quantities are again averaged in 
the same window as above. The variance, albeit clearly decreasing with A^, is in fact still 
so significant as to make an accurate determination of the scaling difficult. Averaging 
over larger time windows the curves become smoother, but such differences persist if we 
use a time window which is small compared to the time scale of the relaxation itself. 
In short the intrinsic finite A^ fluctuations from realization to realization in the (A'^ 
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dependent) relaxation time are still so large for of this order as to limit significantly 
the determination of the average behaviour from a single realization. 

We thus consider a simple ensemble average, over realizations of the initial 
conditions. While we could combine time averaging and such an ensemble average, 
we choose not to do so as this may complicate the interpretation of our result. More 
precisely, if we perform a time average, we would need to check carefully for any possible 
dependence of our results on the chosen averaging window. We will explore below 
in some detail the relation between time averages and ensemble averages over initial 
conditions. 
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Figure 10. Evolution as a function of time (left panel) and (/)22 (right panel), 
averaged in all cases over an ensemble of simulations starting from initial conditions 
which are realizations of cold uniform initial conditions (i.e. i?o = 0). The number of 
particles N and number of realizations averaged over in each case is indicated in the 
panel. The error bars shown are derived (see text) by determining the same quantities 
in two randomly constituted sub-ensembles. 



Shown in Fig. [10] are plots of 0ii and 022 averaged over the indicated number of 
realizations (and without any time average) for each of the indicated particle numbers, 
for i?o = 0. The error bars in this plot have been estimated by dividing randomly the 
realizations into two subsamples and recomputing the average in each of them (i.e. the 
error bar corresponds to the difference in the two averages). 

Using these results we now determine the scaling with A^. Shown in Fig. [TT] is a 
plot of of which treiaxi the characteristic time scale for relaxation, as a function of 
estimated from each of the curves for and 022- We have determined the value of 
treiax in each case as that at which the order parameter reaches half its "plateau" value 
(i.e. in the QSS), i.e., we estimate the value of the parameter which corresponds to the 
approximate plateau and then determine the time at which half this value is attained. 
The error bars correspond to those estimated from those given in the previous figure. 
Shown also are linear behaviours, which in both cases provide a good fit to the results. 

Shown in Fig. [12] are the ensemble averaged evolution of 0ii for the three other 
initial conditions, for the same four values of N . The determinations of the relaxation 
times, for each case where this is possible by the same method as used above, are shown 
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Figure 11. Relaxation time as a function of N , estimated as described in the text 
from the ensemble averaged evolutions of 4>\\ and ^22 shown in the previous figure. 
Linear best fit lines are also shown. The error bars indicated are derived from those 
in the preceding figure. 
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Figure 12. Evolution as a function of time of 01 1 averaged over an ensemble of 
simulations started from initial conditions with i?o = 0.1 (top panel), with i?o = 0.5 
(bottom-left panel) and with Rq = \ (bottom-right panel) . The number of particles 
N and realizations in each case is indicated in the panels. 
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Table 1. Estimated relaxation time treiax for different initial conditions, "-"indicates 
that our data does not allow us a determination of this time using our chosen criterion. 
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in Table [U As there are so few points we have not performed the same fitting procedure 
(with estimated error bars) as above, but it is clear that the given values are consistent 
with a scaling of the relaxation time linear in A^. In the case R = 1, however, we cannot 
deduce any reliable estimate of the scaling with A^, as we can just see the onset of the 
relaxation for A^ = 100 but not in the other cases. 

This last curve and the data in Table [1] allow us to see more quantitatively the 
dependence of the relaxation time on the initial value of Ro also. At fixed A^ we see 
that, between Rq = and Ro = 0.5 the estimated relaxation time increases by a factor 
of about eight. These considerable differences translate into a very different appearance 
to the curves: in the case of i?o = the "QSS plateau" is much less visible as there is 
only a very small separation between the time scales for the establishment of the QSS 
(~ 10^) and the onset of relaxation. 

The exact definition taken here for the relaxation time is somewhat arbitrary — 
we could equally consider the time at which 0ii deviates by 10% from its plateau value, 
or, say, reaches 10% of this value. Because the relaxation is very slow — to show the 
evolution of (pu we must plot it as a function of the logarithm of time — such definitions 
would give enormously different results for the estimated time (differing by two to three 
orders of magnitude). Equally we see from Fig. [11] that if we use 022 rather than 
employing the same criterion we obtain times differing by an order of magnitude. That 
this factor indeed changes only the overall normalization of the times, and not their 
scaling with A^, is evident from the fact that, as can be seen by eye, the curves in the 
decay phase can be superimposed on one another well by a translation parallel to the 
time axis. 

It is interesting to see if a simple functional behaviour can be fit to the decay of 
the order parameters. Shown in Fig. [13] are best fits to two simple functions for the case 
of initial conditions Rq and A^ = 100, for which we have the best statistics. We have 
restricted to the range t > 10^ to cut out the initial (violent) relaxation phase. One 
employs a hyperbolic tangent given by 

^{1 - tanh[a^(logt - logt.e^a.)]} (23) 

in which, therefore, treiax corresponds to the time estimated above. The best fit values 
of the parameters are ipQss = 0.24, treiax = 10"^'^ and ah = 1-4. The other is a stretched 
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Figure 13. Evolution as a function of time of for t > 10'^, calculated over an 
ensemble of one thousand realizations of initial conditions with Rq — sampled by 
N — 100 particles. A best-fit to a hyperbolic tangent [see Eq. (|23|) ] is shown as a 
dashed line, while the solid line is that for a stretched exponential form [see Eq. ([24]) ]. 
The error bars are the same as those in Fig. [TUl 



exponential form: 

0Q55exp{-[-^]"=} (24) 

relax 

and gives the best-fit values (pQss = 0.26 t^eiax — 10'^''' = 0.56. The second 

function is clearly a significantly better fit. We note that the former function has been 
shown in |T8] to give a good fit to the temporal evolution of the magnetisation during 
relaxation in the HMF model. Stretched exponential relaxation, on the other hand, is 
observed in a range of physical systems, and notably in the relaxation of structural and 
spin glasses (see, e.g., [35]). 

We draw attention to one important feature of these results which introduces a 
systematic uncertainty into them, which could only be reduced by doing significantly 
larger simulations: in principle the intermediate QSS is independent of the number of 
particles A^, i.e., we are estimating the dependence of the relaxation time of a state 
which is, up to fluctuations, A^-independent; in practice it is clear in our data that there 
is some residual N dependence in the QSS at the N we are simulating — the "plateau" 
in the curves of the time evolution of our order parameters do not exactly coincide. As 
we have seen that there is clearly a significant dependence of the lifetime on Rq, which 
we interpret to be one on the intermediate QSS rather than the initial condition itself, 
it is possible that the A^ dependence we measure at fixed Rq is due to, or partially due 
to, this residual A^ dependence of the QSS. We believe, however, that such an effect, 
if present, is probably negligible: the differences in the QSS "plateau" at given Rq for 
different A^ are very small compared to the differences between the QSS over the range 
of Rq, and further, for the larger A^, the QSS do appear to converge. This is even the 
case for Rq = 0, where the A^ dependence in the "plateau" is most evident. In this 
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case, as we mentioned above when we discussed our initial conditions, an intrinsic N 
dependence of the QSS might be anticipated: as — t- oo the evolution becomes singular 
at t = 1, and the evolution at finite is regulated by the fiuctuations about a uniform 
distribution which are dependenlQ- That such intrinsic N dependence is weak, if 
present at all, is also indicated by the absence of visible N dependence of 022 in Fig. [101 

4.4- Relaxation and fluctuations in the QSS 

Analytically the relaxation towards equilibrium of systems with long range interactions 
may be described by kinetic equations, derived for example from the BBGKY hierarchy. 
In practice these equations are intractable, and despite many attempts to develop 
appropriate approximation schemes which might make them tractable, there are 
really no solid results which allow us rigorously to model analytically the detailed 
phenomenology of relaxation observed in numerical simulations, and determine for 
example the observed A^ dependence of the relaxation time. 

Inspection of our results for the temporal evolution of the parameters (pu and 022 
lead to one simple observation: the relaxation time appears to be correlated with the 
amplitude of the fiuctuations about the relevant QSS, i.e., the smaller the fiuctuations 
in the QSS, the longer is its lifetime. While this is somewhat trivial when we consider 
a given initial condition (i.e. Rq) at fixed A^ — in postulating that there is a QSS 
we mean that the fiuctuations about it are A^ dependent (and decaying with A^) — 
it is not evident that this should be so for the different i?o at fixed A^. Theoretically 
such a correlation might not be surprising — in kinetic theory approaches the leading 
corrections to the collisionless (Vlasov) limit are, in perturbative approaches, sourced by 
fiuctuations about the QSS (see, e.g., contributions of P.H. Chavanis, and of F. Bouchet 
and J. Barre in [T7].). 

Such a trend can be seen a little in Fig. [71 although in this case it is greatly obscured 
by the time averaging (i.e. it is much clearer if one plots a single realization in each 
case, which we have not done here). It is shown clearly to be present by the results in 
Figs. [H] and [151 The first shows the standard deviation, 0"^^^ , of 0ii as a function of 
time, estimated in the indicated number of realizations of initial conditions Ro = 0, for 
each of the different values of A^ indicated. The error bars in the plot correspond to the 
spread in cr^^^ when it is estimated in two sub-ensembles defined by randomly dividing 
the realizations into two. As remarked above the fact that a^-^-^ decreases with A^ — and 
thus, given that the lifetimes of the states have been observed to increase with A^, that 
there is a correlation of the lifetime with their amplitude — is not surprising: it simply 
means that the fiuctuations about the QSS are, predominantly, due to finite A^ effects 
which will vanish as A^ — )■ oo. We note that the amplitude of a^^^^ in the approximate 
"plateau" region — corresponding to the QSS - scales as 1/y/N, i.e., as they would if 
the fluctuations of (pu is the sum of A^ uncorrelated contributions from the A^ particles. 

For the analogous 3D problem — evolution from cold uniform initial conditions — the precise N 
dependence of the virialized QSS state has been determined numerically in |34j . 
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Figure 14. Evolution as a function of time of ct^jj, the standard deviation of 
calculated in a set of simulations starting from independent realizations of initial 
conditions with Rq — 0. The different curves correspond, as indicated, to different 
values of N and numbers of realizations. The error bars are derived by randomly 
dividing the set of simulations in each case into two subsets. The amplitude of cr^jj at 
any time clearly decreases as N increases. 
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Figure 15. Evolution as a function of time of cr^m the standard deviation of 
calculated in a set of simulations starting from independent realizations of initial 
conditions with the indicated values of Rq, sampled with N = 100 particles (left 
panel) and N — 800 particles (right panel). The number of realizations used in each 
case is indicated. We observe that in each panel the amplitude of an in the QSS phase 
is apparently correlated with the duration of this phase, larger an being associated 
with a shorter relaxation time. In the upper plot the relaxation to thermal equilibrium 
is reflected in the convergence of an for different Rq at later times. 



Shown in Fig. [15] is the same quantity but now for the different values of Rq, at 
two different fixed (A^ = 100 and N = 800). In both cases we see clearly (except 
perhaps for the lowest curve in the lower panel, which is noisier due to the much smaller 
number of realizations) that the amplitude of fluctuations decreases as Rq increases, 
i.e., that the amplitude is (inversely) correlated with the lifetimes we have observed for 
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these states. 

We note that these figures giving the behaviour of the variance of our macroscopic 
parameters also contain alot of other useful information beyond the correlation we have 
just observed. Indeed these curves themselves show very clearly the different time-scales 
in the dynamics: the first period of "violent relaxation" is clearly identifiable by a very 
large variance, which decays on a time scale of order several tens of dynamical times; 
this is followed by an approximately stable value depending on the QSS (Fig. [15]), which 
then evolves on a much longer time scale, dependent on and increasing with A^, towards 
a value which is independent of the initial state (i.e. thermal equilibrium). 

These results also allow us to conclude more about the meaning of our quantitative 
results for the scaling of the relaxation time, which have been calculated using the 
ensemble average: the systematic decrease with of the variance of 0ii (and, we 
have verified, of ^22) implies that such an ensemble average, for sufficiently large N, 
can indeed be interpreted consistently to give the macroscopic behaviour of a single 
realization in the ensemble. Thus, as we have been implicitly assuming, we can 
indeed take our determined relaxation times to represent those of single realizations, 
at sufficiently large N. 

It is interesting to go a little further and consider what the relation is, at finite (but 
large) N, between the fiuctuations in the ensemble average and the temporal fiuctuations 
in a single realization. Indeed if, as we have postulated above, there is a real correlation 
between the amplitude of the fiuctuations measured in the ensemble average and the 
lifetime of the corresponding QSS, it must be that these fiuctuations measured in the 
ensemble bear some close relation to the temporal fiuctuations in the same parameters 
in a single realization. That this is the case, to a good first approximation, can be 
seen from Figs. [16] and [TT] which show exactly the same quantities as in the previous 
two figures, but that the standard deviations are calculated in one hundred time slices 
equally spread over a time window of width At = 10 centred on the indicated point. We 
see the same quantitative behaviours as in the previous plots, and even, in particular 
for = 800, quite good qualitative agreement. 

To allow further more detailed comparison, in Fig. [18] and Fig. [19] are shown, for 
Rq = (left panels) and Rq = 1 (right panels), and = 200 in both cases, the 
histogram of the values of 0ii, at the indicated times measured, in Fig. [18], in one 
hundred simulations from realizations of the same initial conditions, and, in Fig. [TS] in 
one hundred snapshots in a window of width At = 10 in a single simulation from one 
realization of the same initial conditions. Comparing the fours panels in the two figures 
one by one, we see that, although the fiuctuations in each case are clearly not sampled 
from an identical distribution, the agreement is strikingly good: not only, as expected 
from what we have already seen above, do the averages and variances agree well in each 
case, but the general shape of the histograms, which is quite different in each QSS, 
resemble one another strongly. The results are also clearly in line with the conclusion 
drawn above for what concerns the relaxation to thermal equilibrium: at t = 10^ we 
see that the cold initial conditions have relaxed to a distribution centered on the value 
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Figure 16. Evolution as a function of time of cr^u , the standard deviation of 4>ii, for 
the same cold initial condition and values of N as in Fig. [TU but calculated now by 
sampling the value of 01 1 at one hundred equally spaced intervals in a temporal window 
of width At = 10 centred at the indicated time. We see that values are comparable to 
those in Fig. [l^l and show the same trend with N. 
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Figure 17. Evolution as a function of time of the (7^^^ , the standard deviation of 0ii, 
for the same cases as Fig. 1151 but calculated now by sampling the value of at one 
hundred equally spaced intervals in a temporal window of width At = 10 centred at 
the indicated time. The results are very consistent with those in Fig. [T5] 



in thermal equilibrium, while for the case Rq = 1 the system is still in a QSS but has 
evolved very slightly towards equilibrium. 

It would be interesting to develop this study with greater statistics, varying the 
width of the time window to see how good agreement can be obtained, but we will not 
do so here. Such a study is related to the fundamental (and so far unanswered) question 
as to whether the properties of QSS may be determined by averaging over an appropriate 
(non-equilibrium) ensemble, determined by the initial conditions. The theory of violent 
relaxation formulated by Lynden-Bell, for example, postulates an answer to this question 
[36]: the appropriate ensemble is that of all configurations corresponding to a phase 
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Figure 18. Histograms of ttie values of (pu measured at ttie indicated time in one 
hundred simulations started from independent realizations of initial conditions with 
the indicated values of Rq, for N = 200 particles. The dashed line indicates (pu = 0, 
the value at thermal equilibrium. For Rq = the relaxation of the system from the 
QSS to thermal equilibrium is clearly visible, with both the central value and shape 
of the distribution evolving. For i?o = 1 we observe, in line with the results above, 
that the system is still in a QSS but that the distribution has started to shift slightly 
towards the equilibrium one. 



space distribution function permitted by the (collisionless) Vlasov dynamics. If this 
theory were correct (which is not the case for this system [37]) we should perform our 
ensemble average over such configurations rather than the more restricted one we have 
considered. 



5. Conclusions and discussion 
5.1. Summary 

Our primary aim in this paper has been to establish and characterize more fully than 
in the previous literature the relaxation to thermodynamic equilibrium of one of the 
simplest toy models for long-range interacting systems: equal mass self-gravitating 
particles in one dimension (or infinite sheets in three dimensions). Compared to the 
much studied HMF model, notably, the basic properties of this model have remained 
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Figure 19. Histograms of tlie values of ipu measured at one liundred equally spaced 
times in a temporal window of width At = 10 centred on the time indicated in each 
panel. The simulations are from the same cases as in Fig. [181 We observe a qualitative 
agreement between the amplitudes and shapes of those in Fig. [181 



somewhat unclear, and indeed marked by some controversy in the literature. The 
novelty of our work compared to previous studies is not just that we do more and 
larger simulations from a broader range of initial conditions, but that we have identified 
a tool which is very useful to characterize the evolution of the system: the measurement 
of appropriately normalized moments of the distribution function which characterize the 
"entanglement" of the one particle distribution function in configuration and velocity 
space. This is particularly appropriate as a measure simply because the thermal 
equilibrium has the property that such entanglement is absent while, we have shown, 
in any other stationary solution of the Vlasov-Poisson equations such entanglement is 
present. We note that this result, which we showed to be valid for any interaction 
in one dimension (but, as noted, excluding periodic systems like the HMF), can be 
generalized easily to three dimensions if we restrict to stationary solutions which have 
radial symmetry in space and velocity. This suggests that these "order parameters" 
may also be useful indicators of relaxation in much more general, and perhaps, as we 
discuss below, even useful tools for understanding the mechanisms of such relaxation. 

With the aid of these macroscopic measures, we have shown in our numerical study, 
of a range of simple "waterbag" and cold initial conditions, that the system manifests the 
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behaviour thought to be generic in long-range systems: there are essentially two phases 
in the evolution with two completely different time-scales. An initially non-stationary 
state evolves first, on timescales characterized by the "dynamical time" tdyn (roughly 
the crossing time of a particle in the mean-field due to the others), to a QSS, an out 
of equilibrium state, which then evolves on a much longer time scale, dependent on the 
number of particles, to thermal equilibrium. In other words it is reasonable to suppose 
that the system is ergodic (and mixing) on these very long time scales, but not so on 
the shorter time scales. Further we can identify clearly that the QSS resulting from 
different initial conditions (i.e. different values of Rq) are very different macroscopically, 
characterized by very different phase space entanglement. 

Focussing on the the N dependence of the relaxation, averaging over a very large 
number of realizations to average out the fluctuations, we have concluded that the 
characteristic time scale for relaxation behaves, to a very good approximation, as 

trelax ^ f QSS N tdyn (25) 

where fqss is a numerical factor which depends on the initial condition, which we have 
denoted in this way as we expect that this dependence is not strictly on the initial 
condition but on the QSS which results from it. We have seen that this prefactor 
increases as Ro does, by about a factor of ten between Rq = and Rq = 0.5, and 
approximately a further factor of ten for Rq = 1. We have noted that the overall 
normalization of Jqss is rather arbitrary, as it depends greatly on the exact criterion 
used to define the relaxation time-scale. Given that the evolution towards zero of ^n, 
which is what we have used to determine this time scale, is in fact well fit by the simple 
functional behaviour as a function of the time on a logarithmic scale, the normalisation 
of fqss can differ by two orders of magnitude by a trivial change in its definition. More 
specifically we have seen, that in the case where we have accumulated the greatest 
statistics allowing us to constrain the temporal evolution, a very good fit to our order 
parameter (pu is obtained to a stretched exponential form. 

Although the relaxation of this system, and in general long-range interacting 
systems, is not well understood, we can say that this finding of a linear scaling — besides 
the fact that it is, as we will discuss below, in line with less complete previous analyses 
— is not a surprising result: such a scaling can be anticipated both on the basis of 
simple naive estimates of the effects of collisionality, as well from general considerations 
based on kinetic theory. 

A simple "derivation" of this scaling, along the lines of those applied originally 
by Chandrasekhar (see [3H] or [39]) to obtain such an estimate for 3D self-gravitating 
systems, may be given as follows. Relaxation is in principle due to the discretisation, 
in N particles, of a continuous mass distribution. Let us suppose that this latter field 
density varies spatially on a scale, i. The typical fractional change in the velocity v 
of a test particle due to its interaction with one (localized) particle, rather than the 
continuous mass distribution, can be estimated as ~ gi/mv"^. As it crosses the system 
(in a time ~ t^yn) such a particle will interact with all particles. Assuming the 
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contribution from each particle adds incoherently, one estimates 




(26) 



for the normalized variance of the velocity in t^yn- Scaling with at fixed mass and 
energy (and fixed i) requires g/m ~ and therefore ^ ~ A^~^. It follows that the 
relaxation time scales linearly with in units of the tdyn- A slightly different argument 
leading to the same result may be found in |6], and a more developed analysis in |10]. 
In the framework of approaches based on kinetic theory, a linear scaling is obtained as 
coUisional terms arise as the leading corrections in a formal expansion in 1/iV which 
gives the collisionless (Vlasov) limit at leading order (see, e.g. [ITl H2| H3]). 

This scaling linear in is to be contrasted with the case of the scaling observed for 
the life-time in QSS in the HMF, proportional to A^^'^. While the exponent found is not 
understood, the fact that it is larger than unity is consistent with these considerations 
as this result applies for spatially homogeneous QSS (which are possible in the HMF 
because of its periodicity) for which it has been shown that the leading collisional term 
vanishes (see, e.g., contributions of P.H. Chavanis, and of F. Bouchet and J. Barre in 



We note that our study suggests also that the "order parameters" we have defined 
and studied may be relevant quantities for understanding relaxation in this and other 
long-range systems. Indeed in all cases we have observed that, at sufficiently large A^, 
these parameters start from a non-zero value in the initial QSS and evolve monotonically 
towards zero, i.e., the relaxation of the QSS can apparently be described as a process of 
progressive "disentanglement" of the one particle phase space density. In this respect the 
very different, much less efficient, relaxation observed in the HMF might be interpreted 
as a result of the absence of such entanglement in spatially uniform QSS. Further, in 
the case where we have enough statistics to provide a precise fit to the evolution of the 
parameter (pu, we found that it is well fit by a simple stretched exponential form. It 
would be interesting to see in further study whether this fit is more than an approximate 
numerical fit for the case we have studied, and, if so, whether the exponent characterizing 
it is the same or not. As we have remarked such a functional form has been observed in 
other slowly relaxing (e.g. glassy) systems and theoretical tools derived in this context 
to understand relaxation may be relevant. In [35], for example, this behaviour is linked 
to the existence of a fractal structure in a bounded accessible region of phase space. 

5.2. Comparison with previous literature 

Let us now turn finally to compare our findings in greater detail with those in the 
previous literature. 

An early numerical study by Hohl and Broaddus \44\ which concluded a relaxation 
time proportional to N'^tdyn was found to be incorrect by two groups, who studied the 
problem in greater detail (and with greater numbers of particles). However, these groups 
found confiicting results: Miller et al. found no evidence for relaxation at all to thermal 



Relaxation to thermal equilibrium in the self- gravitating sheet model 



29 



equiibrium in their simulations while Luwel et al. [25] found relaxation on a time 
scales even shorter than Ntdyn- Further study (see [5l[27], which also contain a detailed 
synthesis of the previous literature) by Miller et al. concluded that the discrepancy 
was related to the very specific initial condition studied by the other group. Studying 
this case in detail they found that it indeed appears to thermalize very rapidly, but 
some further, but not completely conclusive analysis of the evolution at longer times, 
suggested that this thermalization was not complete. 
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Figure 20. A "counterstreamed" waterbag initial condition in phase space with 
i?o — 0.3, sampled with N — 400 particles. 
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Figure 21. Density profile (left panel) and velocity distribution (right panel) obtained 
at i = 100 from a counter-streamed initial condition with TV = 100. An average 
over 250 simulations from independent realizations of the initial conditions has been 
performed. The solid lines correspond to the values in thermal equilibrium, Eq. (U). 



To determine whether these cases are consistent with our findings — and see 
whether our analysis using the parameters 0ii and 022 can throw light on these 
previous findings — we have resimulated the relevant initial conditions. These are 
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"counterstreamed" waterbag initial conditions, an example of which is shown in Fig. 
We have simulated a range of such initial conditions, in particular the cases (one of 
which is that shown in the figure) considered by [2^ and [5]. Shown in Fig. [2T]are the 
density profile and velocity distribution at t = 10^ obtained starting from a realization 
of initial conditions like those shown in Fig. [201 but with N = 100. We see that the 
profiles indeed agree very well with the equilibrium ones. In Fig. [22] is shown the 
evolution of (pn as a function of time for the indicated values of N averaged in each 
case over the number of realizations indicated. We observe that, although small and 
fluctuating, its value is clearly on average non-zero, indicating that the state, despite the 
good agreement of the density and velocity profiles, is not in fact an equilibrium. Just 
as in the cases we studied we see clearly the relaxation towards equilibrium at longer 
times, and indeed that the characteristic time increases on A^. Although we haven't 
done the more extensive study required to determine precisely this dependence, the 
results are quite consistent with Eq. f[25|) with a value of /qss of order that found for 
the case Rq = 0. 




t 



Figure 22. Evolution as a function of time of 0ii from a counterstreamed waterbag 
initial condition, averaged over the number of realizations of the initial conditions and 
particle numbers indicated. Despite the indications of the previous figure, we observe 
clearly that relaxation to thermal equiibrium has not taken place at t — 100. 

This case illustrates the usefulness of the parameters (pu and 022 as discriminants of 
relaxation: indeed we have just seen that the single measure of (pn is sufficient to discard 
the interpretation of Luwel et al. [26] of their results. This is simply because they are 
physically very appropriate indicators, for the reasons we have explained in introducing 
them: the property they probe — of entanglement of the phase space distribution — 
is one which must evolve significantly during relaxation, because the phase distribution 
must become separable. While (pu and 022 being zero does not imply thermalization, 
of course, we have not found a single QSS, despite exploring a broad range of initial 
conditions (considerably more extended that those reported here) in which they are both 
zero (within the uncertainty of fiuctuations), i.e., the only states we have found in which 
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they are both zero are states which we have concluded, using a range of other measures, 
are indistinguishable from the equilbrium state of Rybicki. It is not difficult, on the 
other hand, to find initial conditions which lead to a QSS in which (pu ^ or 022 ~ 0. 
Indeed for the waterbag initial conditions we have studied both and 022 actually 
change sign as Rq varies over the range we have considered, and one can thus find by 
trial and error the appropriate Rq which make them zero individually. 
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Figure 23. Histogram /(e) of individual particle energies e measured at the indicated 
times starting from counter-streamed initial conditions sampled with N = 100 
particles. The curves are averaged over 250 realizations. In thermal equilibrium /(e) is 
indistinguishable from the one measured at the latest time shown. The result confirms 
that relaxation in fact occurs on time scales similar to those observed for the simple 
waterbag initial conditions. 

Another evident quantity to measure, which we have in fact considered 
systematically but have not reported in detail, is the distribution /(e) of the individual 
particle energy e. This is in fact generally a better discriminant of relaxation than either 
n{x) and /(f), i.e., we have found that in quite alot of cases n(x) and f{v) are not easy 
to distinguish from the equilibrium profiles, but that /(e) allows one to see more clearly 
that one is indeed not in the equilibrium state. An example is the counter-streamed case 
just considered above. In Figs. [23]are shown, for = 100, the evolution of the ensemble 
averaged /(e) at a few different times. We have not plotted the equiibrium curve, as it is 
indistinguishable from the measured curve at the final time shown. One can see clearly 
see that, despite the good agreement of n{x) and f{v) shown in Fig. [211 the system is 
not in equilibrium at the early times: /(e) has a clearly visible excess of particles at high 
energies compared to that at the much later times at which the evolution of 0ii indicated 
relaxation (and /(e) indeed approaches very accurately its predicted equilibrium form). 
While such a measure of /(e), averaged over a large ensemble of realizations, can, in all 
the cases for which we've studied it, clearly discriminate relaxation, the use of just (pu 



Relaxation to thermal equilibrium in the self- gravitating sheet model 



32 



(and possibly ^22) is an extremely efficient short-cut to "diagnose" relaxation. 

Subsequent to |5j, in the nineties, Tsuchiya et al. reported an analysis of larger 
and more importantly, longer, simulations in order to clarify the issue. A ffist paper [6] 
they reported the evolution of a rectangular waterbag initial condition corresponding to 
our case Ro = 1, and reported a detection of relaxation to thermal equilibrium. These 
authors made a distinction between two time scales of relaxation: one of "microscopic 
relaxation" , the other for "macroscopic relaxation" . These are identified, and both found 
to be proportional to A^, by considering the evolution of the mean standard deviation 
of the particle energies averaged over a time window T from their equipartition value. 
The former is estimated from the slope at short time of this function, and the latter 
from the position of "peaks" which are observed to occur at much longer times. While 
the latter is interpreted in terms of of macroscopic relaxation in the sense we have used 
here, the former is interpreted as a time scale on which particles sample the energy 
distribution but on which there is no macroscopic evolution. The justification for these 
interpretations are not clear, and no direct comparison with the equilibrium distribution 
derived by Rybicki, Eq. @, has been reported which might show their correctness. 
Indeed both a subsequent article by the same authors [8] and a study by Yawn and 
Miller [9] place in doubt the correctness of the interpretation in terms of relaxation. 

Nevertheless, in light of the results we have given here, it would be reasonable to 
infer that the results given by Tsuchiya et al. in [6] are indeed correct, at least for 
what concerns the N dependence of the relaxation. Further comparison could of course 
clarify the relation of the behaviour of their measured quantities and the macroscopic 
relaxation as we have probed it here (and should be much easier for the shorter lived, 
smaller Rq, initial conditions rather than the case Rq = 1 studied by these authors). 
We do not believe, however, that there is any clear basis for either an operational or 
physical distinction between "microscopic" and "macroscopic" relaxation as described 
by these authors: as we have discussed there is an arbitrariness in the definition of 
the relaxation time because of the very slow nature of this relaxation. As we have 
noted, we could easily, for example, have obtained here estimates of the relaxation time 
differing by several of orders in magnitude in their prefactor, just like the two different 
time scales determined by Tsuchiya et al. [6], by using slightly different definitions, or 
choosing to use a different order parameter. This point can be illustrated by considering 
the evolution of /(e) for one of the cases we have considered: shown in Figs. [2l]is this 
quantity for the case R = 0.1 and A^ = 400, averaged over 60 realizations. While we have 
associated (see Table [1] above) the time scale 7 x 10^ to the relaxation in our analysis, 
one can discern by inspection of these figures significant evolution (in particular of the 
initially clear "core-halo" structure) in /(e) already by t = 10^'^, i.e., there is evolution 
of the energy distribution on the time scale of "microscopic" relaxation (of order Ntdyn) 
identified in [6]. While it is possible that there are different time scales associated to 
different physical processes as argued in [6], it seems a more plausible interpretation 
to us to suppose that there is single physical relaxation process leading, albeit very 
slowly, to macroscopic relaxation of the system, and to characterize this relaxation by 
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a function and the scaling of its parameters with A^. In this respect it is interesting to 
note that the specific stretched exponential form we fitted to the temporal behaviour 
has the known property [16] that it is can be written as a weighted integral over simple 
exponentials (i.e. it can be interpreted as arising from the superposition of an infinite 
number of relaxation processes each with a single characteristic time). 
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Figure 24. Histogram /(e) of individual particle energies e at the indicated time, and 
averaged over 60 simulations from realizations of simple waterbag initial conditions 
with Rq — 0.1 and N — 400. The curve at the latest time coincides well with that 
expected in thermal equilibrium. The onset of relaxation is already visible at a time 
of order 10^, almost two orders of magnitude smaller than the time determined in 
Table [H 



In [8] Tsuchiya et al. have also described chaotic "itinerant" behaviour of 
these systems, starting from the same {Rq = 1) initial conditions i.e., in which the 
system shows apparently stochastic macroscopic behaviour. In our analysis this would 
correspond to such behaviour for the parameters (pu or 022- While we have seen that 
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there are indeed very significant fluctuations in these parameters, which correspond to 
very signiflcant differences in the "macroscopic" evolution of these systems, we have 
studied carefully their dependence on and found them to decay monotonically. The 
results of [8J were obtained for = 64, a range in which we still see fluctuations of 
011 which are order unity. Only when we reach of order several hundred do we 
see these fluctuations diminish signiflcantly so that the macroscopic trajectory of the 
system becomes quite localized. We thus beheve that as N increases these effects will 
becomes negligible, even on the time scales on which relaxation occurs, and an effectively 
deterministic macroscopic evolution will occur. 

It is interesting to compare our results also to those of Yawn and Miller [1], who 
have analyzed in detail relaxation in a version of the sheet model in which there are 
sheets of different masses. In this case the relaxation towards thermal equilbrium may 
be clearly distinguished by testing for equipartition of the kinetic energy, and the 
associated spatial segregation of the mass populations. In simulations starting from 
waterbag type initial conditions with a virial ratio of two, for a range of different mass 
ratios and up to = 128 particles, clear evidence was found in |1| for such relaxation 
using appropriately deflned indicators. Like the order parameters we have employed 
here, these show characteristic behaviours corresponding to the principal phases of 
the dynamical evolution (violent relaxation, QSS phase, relaxation towards thermal 
equilibrium). Although we cannot compare our results directly, we note that the time 
scales observed for relaxation of systems with A^ ~ 10^ particles are quite consistent 
with those we have observed for the equal mass system with initial virial ratio Rq = 1. 
Yawn and Miller [1] also measure temporal correlation properties and flnd weak but 
persisting correlations characterized by a power-law decay (in time), which they interpret 
as evidence for the incompleteness of relaxation. In the present study we have found, 
in contrast, that our principal observables decay in time with a functional form which 
allows the identification of characteristic time scales. Further all deviations of these 
observables from their equilibrium values decrease clearly as A^ increases, and thus we 
have interpreted the associated "incompleteness" of relaxation simply in terms of finite 
A^ effects. It would be interesting certainly to perform a more direct comparison of the 
results in the two models, and in particular to extend the study of Yawn and Miller to 
allow a determination of the A^ dependence of the parameters they study. We note also 
that Yawn and Miller argue that the power-law decay suggests the existence of a fractal 
structure in phase-space, which, as we have been mentioned above, is also proposed 
as an explanation in [35] for the appearance of relaxation characterized by a stretched 
exponential behaviour. 

Let us finally mention some other issues of importance concerning aspects of the 
dynamics of this system which have been treated elsewhere but which we have not 
discussed here. As we have discussed, we interpret our results in line with those of 
many previous studies of this and other long-range systems: the evolution from an 
arbitrary out of equilibrium initial condition is characterized a first phase of relaxation 
to a QSS, interpreted as a finite particle sampling of a stationary solution of the Vlasov 
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equation, on a time scale independent of N, followed by a slow relaxation to thermal 
equilibrium on an A^-dependent time scale. Studies of the single mass sheet model for 
other specific initial conditions suggest that this simple scheme may be too limiting, for 
this model (and possibly, for all such models). On the one hand Reidl and Miller have 
reported numerical results [17] for specific "two cluster" initial conditions which show 
a dependence on N in the time scale for relaxation to a QSS. On the other hand, as 
mentioned in the introduction, Rouet et al. [281 129| have shown, using both particle 
simulations and simulations of the Vlasov equation, for yet other initial conditions 
that "holes" which rotate in phase space may be present after violent relaxation and 
persist on very long time scales. Although it is not evident that there is necessarily a 
relation between either finding and the mechanism of relaxation to thermal equilibrium, 
a study incorporating such initial conditions would certainly be more complete that that 
reported here. Extension of the study reported here to larger still would likewise be 
desirable, despite the extremely rapidly growing numerical cost of such simulations with 
N. 
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